lemur.data <- read.csv('Lemurdata_glm_2.csv')
lemur.data <- read.csv('Lemurdata_glm.csv')
lemur.data <- read.csv('Lemurdata_glm.csv')
# access the variables in the dataframe without having to call the dataframe
# i.e. can print 'height' instead of 'lemur.data$height'
attach(lemur.data)
m2 <- glm(GIparasites ~ sex+year, family=poisson, data=new.dat)
# let's first ask how do the number of GI parasites change over a given year? We probably want to exclude the
#intervention year (2023)
new.dat <- subset(lemur.data,year!='2022') #you should now have only 300 observations
m2 <- glm(GIparasites ~ sex+year, family=poisson, data=new.dat)
summary(m2)
m2 <- glm(GIparasites ~ sex+ (1|year), family=poisson, data=new.dat)
summary(m2)
View(new.dat)
# now that we see the data, let's use a glmm to ask whether the sites are different, but I want to control for year
# and individual ID!
m4 <- glmer.nb(GIparasites ~ location + sex + (1|year) + (1|ID), data=new.dat)
library(lme4) #generalized linear model package used on Monday
library(lmerTest)
library(ggplot2) #ggplot visualization package. We will alternate between this and the base R plot() functions
# now that we see the data, let's use a glmm to ask whether the sites are different, but I want to control for year
# and individual ID!
m4 <- glmer.nb(GIparasites ~ location + sex + (1|year) + (1|ID), data=new.dat)
summary(m4)
# okay sex is significant and since we know it it important, but do not necessarily care about it's effect,
# let's add it as random
m5 <- glmer.nb(GIparasites ~ location + (1|sex) + (1|year)+ (1|ID), data=new.dat)
summary(m5)
# okay sex is significant and since we know it it important, but do not necessarily care about it's effect,
# let's add it as random
m5 <- glmer.nb(GIparasites ~ location + (1|year)+ (1|ID), data=new.dat)
summary(m5)
# okay sex is significant and since we know it it important, but do not necessarily care about it's effect,
# let's add it as random
m5 <- glmer.nb(GIparasites ~ location*sex + (1|year)+ (1|ID), data=new.dat)
summary(m5)
# But if your groups are different, which ones?
# in this instance, we might want to run a post-hoc test. We need to load a new package to look at post-hoc tests
#library(multcomp)
summary(glht(m5, mcp(location="Tukey")))
# But if your groups are different, which ones?
# in this instance, we might want to run a post-hoc test. We need to load a new package to look at post-hoc tests
library(multcomp)
summary(glht(m5, mcp(location="Tukey")))
#let's plot this just for fun!
ggplot(data=new.dat, aes(x=location, y=GIparasites, fill=sex))+
geom_boxplot()+ theme_classic() + xlab("\nLocation") + ylab("Number of GI parasites\n")
summary(glht(m5, mcp(location*sex="Tukey")))
# So let's run a model
m6 <- glmer.nb(GIparasites ~ treatment + (1|sex) + (1|ID), data=lemur.data)
summary(m6)
library(lme4) #generalized linear model package used on Monday
library(lmerTest)
library(ggplot2) #ggplot visualization package. We will alternate between this and the base R plot() functions
lemur.data <- read.csv('Lemurdata_glm.csv')
# let's first ask how do the number of GI parasites change over a given year? We probably want to exclude the
#intervention year (2023)
new.dat <- subset(lemur.data,year!='2022') #you should now have only 300 observations
#-------------------------------------#
#   4. Generalized additive models: GAMs ------
#--------------------------------------#
lemur.data$year <- as.numeric(lemur.data$year)
library(mgcv)
gam1 <- gam(GIparasites ~ s(as.numeric(year), k=-1,bs="tp"),  data=lemur.data, family="gaussian")
gam1 <- gam(GIparasites ~ s(as.numeric(year), k=4,bs="tp"),  data=lemur.data, family="gaussian")
gam.check(gam1)
gam1 <- gam(GIparasites ~ s(as.numeric(year), k=7,bs="tp"),  data=lemur.data, family="gaussian")
gam1 <- gam(GIparasites ~ s(as.numeric(year), k=5,bs="tp"),  data=lemur.data, family="gaussian")
gam1 <- gam(GIparasites ~ s(as.numeric(year), k=2,bs="tp"),  data=lemur.data, family="gaussian")
gam1 <- gam(GIparasites ~ s(as.numeric(year), k=3,bs="tp"),  data=lemur.data, family="gaussian")
gam.check(gam1)
gam1 <- gam(GIparasites ~ s(as.numeric(year),bs="tp"),  data=lemur.data, family="gaussian")
gam1 <- gam(GIparasites ~ s(as.numeric(year), k=4,bs="tp"),  data=lemur.data, family="gaussian")
pred.df <- dplyr::select(lemur.data, year)
summary(gam1)
